Emergence of coherent backscattering from sparse and finite disordered media

Coherent backscattering (CBS) arises from complex interactions of a coherent beam with randomly positioned particles, which has been typically studied in media with large numbers of scatterers and high opacity. We develop a first-principles scattering model for scalar waves to study the CBS cone formation in finite-sized and sparse random media with specific geometries. The current study provides insights into the effects of density, volume size, and other relevant parameters on the angular characteristics of the CBS cone emerging from sparse and bounded random media for various types of illumination, with results consistent with well-known CBS studies which are typically based on samples with much larger number of scatterers and higher opacity. The enhancements are observed in scattering medium with dimensions between 10× and 40× wavelength and the number of particles as few as 370. This work also highlights some of the potentials and limitations of employing the CBS phenomenon to characterize disordered configurations. The method developed here provides a foundation for studies of complex electromagnetic fields beyond simple incident classical beams in randomized geometries, including structured wavefronts in illumination and quantized fields for investigating the effects of the quantum nature of light in multiple scattering, with no further numerical complications.

cone 1 . Furthermore, numerical simulations such as Monte Carlo and related approaches 18,[28][29][30] , and T-Matrix methods 31,32 have also been used to model the multiple scattering processes in disordered media.
Besides isotropic samples, the nature of light scattering and transport in disordered samples with structural anisotropy have been the subject of multiple numerical [33][34][35] , theoretical 36,37 and experimental 34,35,37,38 investigations. Many non-biological samples such as nematic liquid crystals, porous fibers, and porous semiconductors as well as biological samples such as muscle, bone, and skin tissues may exhibit such structural anisotropy, causing them to show non-symmetric near-and far-field scattering patterns or anisotropic transport of light. Transport of light in anisotropic samples is found to be adequately captured via an anisotropic diffusion model 36 for cases where the propagation length is several times larger than the transport mean free path. Consequently, characterizing such diffuse transport can determine the diffusion constant, the mean free path tensor, and the extrapolation length 38 . In addition, the multiple scattering interference phenomenon in such systems can produce angularly anisotropic behaviors in weak localization of light, i.e., anisotropic enhanced backscattering 39,40 , which can be effectively captured considering enhanced backscattering models with generalized anisotropic mean free paths 39 . In this regard the CBS phenomenon has been used as a static and accessible technique to characterize the optical anisotropy of biological samples such as beetle scales 41 .
In this paper, we use a microscopic numerical scalar wave analysis to examine the formation of CBS peaks in scattering from finite-sized, wavelength-scale random media. The inherent simplicity of the calculation method allows us to systematically study how different parameters of the scattering object can affect the line shape of the backscattered enhancement, and specifically to understand the limits on the formation of the backscattering cone for small numbers of scatterers and sparse samples. Additionally, this approach allows a straightforward way to calculate the effects of density factor, sample size, absolute number of particles, and iteration factor. Here, we first consider scattering from small scatterers positioned within a rectangular space with different edge sizes ranging from 10 to 40 free-space wavelengths, with excitation in the form of a monochromatic plane wave. This approach allows us to easily extend the computation to the scattering of arbitrary input states of light, such finite-size excitation beams or other tailored excitations 14,42 . In addition to the emergence of CBS, we analyze the ensemble-averaged specular reflection in the scattering from random samples, a feature typically overlooked in previous numerical studies. We then develop a method to eliminate this coherent specular contribution, and hence enable a clean observation of the pure CBS peak under normally incident illumination.
While a scalar analysis does not capture polarization effects [43][44][45] (e.g., the polarization opposition effect 43 ), it allows us to study more populated scattering volumes. Here, we analyze volumes consisting of up to 20,000 particles, much higher than previous studies which computed exact solutions of Maxwell's equations 46 .

Formalism and modeling
To model the CBS formation, a random medium is typically approximated by a half space in order to simplify the calculation 1,3,47 . This approximation is accurate for cases of centimeter-sized samples under external optical illumination, as the mean free path length of the wave inside such disordered media is much smaller than the physical size of the scattering sample. However, such models do not provide insights into the behavior and physics of wave scattering in the case of smaller random configurations, i.e., those with dimensions in the order of multiple wavelengths or comparable to the mean free path length, as we study in this work. Figure 1 depicts a sketch of the finite-sized random media examined in this study together with the direction of illumination with respect to the chosen coordinate system. A total of N p isotropic small particles (scattering centers) are assembled inside an imaginary rectangular domain, where the center of each particle is randomly Figure 1. Simulation setup composed of N p particles randomly arranged inside a rectangular region with dimensions L x , L y , and L z . The structure is illuminated by a plane wave propagating in the xz-plane with the incident angle of θ i . The area between the dielectric particles is filled with air and the scattering is calculated in the upper half plane for −π 2 ≤ θ s ≤ π 2. www.nature.com/scientificreports/ positioned. We use uniform statistical distributions to assign all three components of the particles' positions (i.e., x, y, and z); and the process is repeated N r times to realize N r independent random distributions. In experiments, ensemble averaging of the scattered intensity, which results in a significant reduction of granular behavior (speckle), is typically achieved through rotating or moving the sample. This step is executed here, however, by computing fully independent realizations of the composition followed by ensemble averaging over the calculated results. The uniform distribution of the particles' positions (which could be altered to any other distribution) is the only assumption made about the scattering medium, and we do not impose any assumptions on the statistical behavior of the response or the scattering matrix elements. It is noteworthy that since N p specifies the number of equations that needs to be solved at each iteration, changing the distribution of particles, their composition, or the domain size, does not affect the computation load and processing time. We have found that in this model, a typical personal computer with 32 GB RAM is capable of handling N p values in the order of up to about 10 4 ; no cloud computing or advanced computing techniques were required to obtain the results presented here. The structure is illuminated with a scalar electromagnetic plane wave propagating in the xz-plane, incident upon the structure at the angle θ i , as shown in Fig. 1. The far field scattering pattern from the particle cluster for all outgoing wavevectors (depicted as θ s ) is then calculated. For −π 2 ≤ θ s ≤ π 2 the entire upper half-plane region is covered, in which θ s = −θ i corresponds to the backscattering direction and θ s = θ i corresponds to the specular reflection of the incident wave. With the finite size of the sample, the random object also creates scattering in the forward direction (i.e., |θ s | ≥ π 2 ), which is not of interest in our studies. In the back half-plane, as will be discussed below, the observed coherent effects in multiple scattering, which survive the ensemble averaging, are mainly at the specular ( θ s = θ i ) and the backscattering ( θ s = −θ i ) directions.
To calculate the scattered field, we solve for the total scalar field amplitude U(r) satisfying the general Lippmann-Schwinger integral equation given in Eq. (1) and independently solved for each realization. The total field amplitude at r can be constructed from the incident field amplitude U i (r) and the contributions from all scatterers, where G 0 r, r ′ is the background Green's function 48 and k 0 is the free space wavenumber The susceptibility parameter η(r) is defined such that it is locally related to the relative permittivity by ε r (r) = 1 + 4πη(r) 48 . For a collection of discrete point-like scatterers (i.e., delta function approximation), η(r) may be described as a summation over the effective polarizability of particles defined at the center of each particle The free space Green's function G 0 r, r ′ = e ik 0| r−r ′ | r − r ′ is utilized, following the e −iωt convention, with proper adjustments for the regularized Green's function at r = r ′50 . In the final step, the scattering amplitude of the outgoing spherical wave, in the far field domain, is extracted and studied (See supplemental information for further details on solving the Lippmann-Schwinger equation and calculating the scattering amplitudes). We model N p identical dielectric particles (i.e., simulating commonly used SiO 2 particles with refractive index n = 1.5 ) and subwavelength radius of r = 0 2π . Clearly, the most accurate solution may be attained by using the complete Mie scattering theory for a composition of spheres while solving for the fields using the exact Maxwell's equations 32,51 ). The computation time and complexity of such analysis, however, drastically increases with increasing the number of the scatterers and the size of the simulation domain. On the other hand, the scalar approximation 1,6,17,47 enables us to capture many important underlying physical properties of the system, while avoiding unrealistically long simulation times.
Several physical factors must be considered in the choice of the material and the size of the scattering bodies. In particular, if the scattering elements are very far off resonance (i.e., with extinction cross sections much smaller than the physical cross section of the particles), one might expect that millions of particles are required to create a meaningful scattering response, associated with large mean free path lengths 2 .
Along the same line of physical reasoning, we expect to observe meaningful coherent scattering responses and weak localization effects from a smaller number of scatterers when the scattering amplitude is higher. Using Mie theory, the extinction efficiency (i.e., extinction cross section normalized to the physical cross section of the particle) for the dielectric spheres (scattering centers) modeled in this work is found to be 0.22, the amplitude of the electric dipole scattering coefficient is 0.19, and the amplitudes for all other higher order scattering coefficients are smaller than 0.03. We note that the response of each single particle is still well within the off-resonant dipolar approximation. For the lossless configurations examined here the extinction cross section is equal to the scattering cross section and represents the amount of power that the object extracts from the input wave normalized to the incident power intensity.
In the following, we study the formation of the CBS cone from various finite-sized media based on the above formulation. In addition to the conventional parameters applicable in studying semi-infinite random media, such as density factor and angle of incidence, here we also study parameters specifically relevant to finite-sized systems including the absolute number of scatterers and the size of the domain. The ensemble-averaged specular reflection in multiple scattering and a technique to eliminate it is also discussed and the modified results are presented.

Results
We analyze several examples of random geometries (such as that depicted in Fig. 1) following the numerical formulation described in the previous section. In all the examples, the random domain is a cube, i.e., L x = L y = L z = L and density factor is the ratio of the volume of all particles to the volume of the cube, defined as www.nature.com/scientificreports/ in which V S is the volume of each dielectric sphere. The ρ parameter is an indication of how closely the particles are packed. The scatterers are assumed to consist of a low-index material ( n = 1.5 ) with a radius of r = 0 2π . As discussed above, the Mie scattering analysis of individual particles shows that their response lies well within the single-mode dipolar scattering regime and higher order scattering modes are negligible. For all the results provided in Figs. 2, 3, 4, 5, 6 and 7, N r distinct random configurations are simulated, and the results are averaged to reveal the coherent contributions vs. the incoherent background signals. We first assess the effect of disorder averaging. For samples with a very large number of particles (the case typical of most laboratory experiments), usually tens to hundreds of measurements are sufficient for the speckle to average to a diffuse background allowing the CBS cone to emerge 52 . This relatively low number of realizations is not sufficient to observe the weak localization effects in small, disordered geometries. By increasing the number of realizations, however, the localization effects appear for these geometries as well. This important conclusion is carefully examined here. Keeping all the other parameters unchanged, Fig. 2a illustrates the evolution of the averaged scattering profile when N r is increased from 20 to 10,000 in a random sample with N p = 10, 000 , L = 20 0 , θ i = 40 (deg) ≈ 0.7 (rad) . This corresponds to a density factor of ρ = 0.021 or approximately 2% volume fraction. Here, 0 indicates the free space wavelength.
First, we notice that ensemble averaging over 2000 (or even smaller) number of simulations is sufficient to reveal the coherent scattering response in these geometries. The CBS cone is still observable for averaging over 200 samples, however as N r increases the response across all angles becomes less noisy. Another important observation is that the ensemble-averaged specular coherent reflection is present at θ s = θ i in averaging over 200 realizations or more This effect is well separated from the CBS contribution under oblique incidence and no specular reflection is present in scattering from any single realization (Fig. 2b) or small values of N r . Note that the background material here is air everywhere; thus, all the observed specular effects are the direct consequence of embedded coherence effects in multiple scattering, outliving the ensemble averaging step. The CBS ratio is slightly smaller than the theoretical value of two and does not change noticeably by increasing the number of realizations. The finite size of the scattering medium plays a fundamental role in establishing the width of the CBS cone. For finite-sized samples, such as the examined cases here, the probability of having longer scattering trajectories for photons severely diminishes, and thus the width of the CBS cone is anticipated to increase [53][54][55] . This is consistent with our observations here (see Figs. 2, 3, 4, 5, 6, 7) where the finiteness of the samples results in broader CBS cones.
Performing a statistical analysis over multiple implementations of the random medium provides insight into the coherent wave effects inside the random structure; at the same time, it is also insightful to investigate the scattering response for individual realizations. Here, to emphasize the effect of ensemble averaging, the angular distribution of scattering intensity for two sample realizations (i.e., N r = 1 ) is shown in Fig. 2b. These curves demonstrate the speckle patterns in the far field domain, originating from interference of waves due to multiple scattering paths.
Next, we investigate the effect of the angle of incidence on the emergence of weak localization and the degree of enhancement as well as the evolution of the specular fraction of the coherent scattering. A particularly interesting case is when θ i = 0 (deg) , in which the backscattering and specular directions overlap. This is shown in Fig. 3 where we examine four incident angles, θ i = (60, 40, 20, 0) (deg) covering from near-grazing angles to exact normal incidence. Physical properties of the sample are unchanged compared to the results in Fig. 2 (   Fig. 3, where an artificial random surface is added on top of the medium shown in Fig. 1. The angle of incidence is set at θ i = (60, 40, 20, 0) (deg) (different colors). The darker colors corresponding to higher incidence angles. N r = 2, 000 and the maximum value (CBS peak) of all curves are normalized to two. Each curve is moved one unit up compared to the previous one for better visualization.  Fig. 3, the width of the CBS cone and the enhancement ratio is approximately independent of the angle of the incidence except for the case of normal incidence. The fact that the width of the peak only depends on the properties of the sample is consistent with previous studies indicating a direct relation between the width of the CBS cone and the mean free path length inside the disordered medium 1,14,56 . We also note that the width of the cone may be slightly overestimated for smaller values of N r . Interestingly, for the case of normal incidence, we notice that both the line shape and width of the CBS cones deviate from the results for oblique incidence. This is due to the superimposition of the specular reflection and the coherent backscattering at this angle, which makes it difficult to accurately estimate the line shape and width of the CBS cones.
Indeed, an interesting observation in all the numerical results presented so far is the occurrence of specular reflection arising in scattering from the random samples. Experimentally, to avoid the specular reflection, random samples are typically tilted off axis and ensemble averaging over these angles are not of interest 6 . However, here, we closely examine the evolution of averaged intensity over these angles and devise a technique to explicitly eliminate averaged specular peaks while maintaining the localization contributions unchanged. Numerical results in Fig. 3 also indicate that quite contrary to the CBS, specular intensities are very sensitive to the angle of incidence. For near-grazing angles, both the intensity and width of the specular reflection are more pronounced compared to normal illumination.
The emergence of the specular reflection from statistical averaging can be interpreted as the direct consequence of the "artificial interface" present at z = L z . This artificial interface is built by restricting the z position of all particles to be inside the cube. For finite-sized samples, as is the case here, we observe comparable specular  www.nature.com/scientificreports/ and CBS intensities. To eliminate the specular peak, the artificial interface (on top of the block) is replaced by an effective "random surface" stretching over a finite thickness around z = L z . Note that this "imaginary" interface, in planar or random form, has a real existence with respect to the average field which is often considered in theoretical studies of the system 1 . The local position of the artificial random surface at each point in space is determined using semi-Fourier series functions constructed using trigonometric functions. As the result, the new random medium is uneven at the top interface and the coherent properties at the specular angle are expected to diminish. This is indeed the case as shown in Fig. 4, for random samples with N p = 10, 000 , L = 20 0 , and N r = 2, 000 examined under various illumination angles. Compared with Fig. 3, the specular wave is clearly suppressed and is totally undetectable. Notably, the CBS cone under the normal incidence is now wider (comparing Figs. 3d and 4) and the observed cone width is consistent with the widths under oblique illumination as well. Also, the degree of enhancement for normal incident wave modifies to a value slightly below 2 in the new simulation, consistent with the ratios observed at oblique incidence, as expected. Therefore, this approach makes it possible to accurately characterize the coherent backscattering effect, even at normal incident angles. We next examine the line shape of the CBS cones observed for the finite-sized, wavelength-scale geometries, and compare these intensity profiles with theoretical predictions traditionally made assuming large configurations. In the case of a semi-infinite scattering medium, the backscattering cone follows a triangular shape (i.e., linear behavior near the backscattering direction). When the medium is of finite size or finite thickness, the triangular line shape becomes rounded due to the omission of the higher-order multiple scatterings 1,53 . The theory of the CBS from finite slabs was first proposed in using a scalar wave theory 6,53 and has been also studied in more recent works 41 . Generally, the line shape of the backscattering cone depends on the wavelength and the mean free path length inside the random media 1,2,5 . In view of the fact that the CBS enhancement effect is present in the case of finite sparse random geometries with the intensity profile behavior similar to semi-infinite configurations, it is important to examine this similarity more closely. Specifically, we study if the predicted enhancement profile for large configurations can still adequately captures the enhancement line shape for small sparse configurations. To this end, we compare the profile of the localization here with the theoretical model devised by Akkermans and Montambaux 2 with the correction for finite-sized random domains Here, b = L l e , where L is the thickness of the sample and l e is the mean free path length used here as a fitting parameter. k ⊥ is the amplitude of the transverse projection of the vector k i + k s on the xy plane, in which k i and k s are incident and scattering wavevectors, respectively 2 . Figure 5 compares the numerically calculated line shape with a theoretical fit using Eq. (3) for a sample with 10,000 particles and box size of L = 20 0 . The multiple scattering is happening here in the weakly disordered regime above the Ioffe-Regel criterion 57 for the appearance of Anderson localization. We note that the shape of the backscattering cone for this finite, wavelength-scale geometry using the scalar wave analysis can still be closely captured by those analytical predictions (see supplemental information for further discussions). www.nature.com/scientificreports/ To better understand the effects of sample size, number of particles, and density factor on the shape of the cone, in the following we investigate two specific cases where we fix the density factor while changing the box size (Fig. 6), and number of particles (Fig. 7). Figure 6 illustrates the combined effects of density factor and sample size on the characteristics of the CBS cone. The number of particles is set at a reasonably high value of 10,000, while the box size is gradually increased from 10 0 to 40 0 . This is associated with decreasing the density factor from approximately 17 to 0.26% and we conduct the simulations for 1000 (Fig. 6a) and 10,000 (Fig. 6b) independent realizations. Figure 7 shows the combined effects of the absolute number of particles (i.e., N p ) and sample size, for two fixed density factors. In Fig. 7a, density factor is fixed at 0.63% and the box size is gradually increased from L = 10 0 to L = 34.3 0 , corresponding to samples with very few particles (i.e., N p = 370 ) to very large number of particle numbers (i.e., N p = 15, 000 ), respectively. In Fig. 7b, the density factor is set at 2% and the box size is changed from L = 10 0 to L = 25.2 0 .
Several interesting observations can be made from the results presented in Figs. 6 and 7. First, samples with lower density factor are more sensitive to ensemble averaging, so a higher number of independent realizations is required to retrieve the embedded coherent effects. In addition, the width of CBS cone is slightly overestimated for smaller N r , also more evident for sparse samples. Second, the intensity of the specular reflection is directly related to the density factor, as expected. At the studied size range (between L = 10 0 to L = 34.3 0 ) the width of the specular reflection is slightly dependent on the size of the sample, and we observe narrower specular reflections for larger boxes. This width, however, is significantly decreased for sparse samples (Fig. 6). Third, we notice that the coherent contributions in the backscattering, i.e., both CBS and specular reflection, still emerge for even as few as a few hundred particles (Fig. 7a), corresponding to very sparse random samples with densities as small as ρ = 0.0063 . To capture such coherent features, ensemble averaging on many different samples may be required. Fourth, the width of the CBS cone is strongly proportional to the density factor with a slight inverse dependence on the absolute number of particles for these samples. This is consistent with previous studies on the role of mean free path length on the shape of the cone, FWHM ∝ 1 l e 14,56 , indicating longer path lengths for more sparse samples as long as the sample size is much larger than l e . Fifth, the enhancement ratio is closer to the theoretical limit of two for samples with higher density factor, consistent with previous analysis of electromagnetically large random samples 5,27,54 . In the case of constant density factor, increasing the sample volume (i.e., more particles) also significantly increases the enhancement ratio (Fig. 7).
We note that our approach to study the emergence of the CBS cone in small random samples may be utilized to investigate other forms of scalar wave interaction with disordered structures (besides CBS). As an example, illumination with structured profiles, such as those generated using spatial light modulators, can be implemented with no additional numerical complications. We also provide a quantitative study of the behavior of CBS cone emergence, information that is valuable in experimental studies. The number of realizations sufficient to converge the enhancement and the behavior of convergence is not a quantity accessible through techniques relying solely on averaged quantities. Additionally, the fact that the CBS enhancement is present for very small samples and for very low densities suggest the possibility of characterizing such samples using CBS measurements. This technique is already widely in use for characterizing optically thick samples mostly with large opacity, and our results indicate the possibility of extending this technique to study other structures such as thin films (which are among widely used platforms in scientific and industrial communities for various applications). The level of disorder or porosity of such films can degrade the performance of the device or on occasion be utilized to gain a benefit, and thus must be carefully characterized. Studies of such samples with finite illumination and collection apertures is effectively equivalent to studying sparse and small structures (in all three dimensions). Therefore, our studies on CBS cone for finite-sized geometries suggest a useful technique to characterize the spatial inhomogeneity of samples.

Conclusions
In this work we have provided a comprehensive study of the coherent backscattering from finite-sized sparse scattering medium based on a scalar wave analysis. Relying on a simplified form of Maxwell's equations in a multi-scattering medium, we investigated the effect of density factor, number of particles, and size of the scattering medium on different characteristics of the CBS line shape. The density factor was shown to be the prominent factor in determining the width of CBS cone, while both the density factor and the number of particles influence the degree of enhancement in the backscattering direction. The angular distribution and the enhancement factor of the scattering intensity hold valuable information to characterize random samples. The cone shapes are found to be in very good agreement with previous theoretical calculations. Finally, we have devised a technique to exclusively eliminate the specular reflection component from the ensemble-averaged backscattered intensities, enabling the observation of pure CBS line shapes even for normal illuminations of the structure.
While the results of this work are entirely consistent with the known properties of CBS in various systems, it is valuable in designing and interpreting experiments on practical, finite-size systems to be able to compute the detailed scattered field and observe the emergence of the main features of the CBS light. A relevant situation is that of coherent backscattering enhancement in single scattering phenomenon. This effect has been recently studied 58 by making connections between the solution of Maxwell's equations in scattering from single particles and an order-of-scattering formalism. Interference effects, analogous to the ones discussed here, between the conjugate pairs of sequences result in enhancements. It has been shown that the angular width of coherent backscattering enhancement-induced peak is significantly dependent on the particle size and is inversely proportional to this parameter. Similar behavior has been observed in our work (see Fig. 6). To make the connection, a finite-size sparse medium captures some of the important behaviors of a single large scattering object, including the behavior of coherent backscattering enhancement in single scattering. Also, the approach and results developed here Scientific Reports | (2022) 12:22256 | https://doi.org/10.1038/s41598-022-25465-y www.nature.com/scientificreports/ may be most useful as a foundation for the study of multiple scattering in random media when the incident field is more complex, such as highly structured classical light fields or fields with various types of quantum correlations. Given that the scalar fields are calculated for each realization and the statistical averaging is performed through direct ensemble averaging, as opposed to finding statistically averaged fields, it is now possible to study such geometries in the context of quantum optics. Specifically, the spatial profile of the scattering can be fully studied classically while the quantum nature of the studies can be captured in the creation and annihilation of the photons and the photon number distributions.

Data availability
Data that supports the findings of this study are available from the corresponding author upon reasonable request.